using System;
using System.Collections.Generic;
using System.Linq;

namespace MT5Trade.Models
{
    /// <summary>
    /// GARCH(1,1)模型实现类
    /// 实现标准的GARCH(1,1)参数估计算法，为波动率预测提供数学基础
    /// 遵循严格的数学模型：σ²ₜ₊₁ = ω + α × rₜ² + β × σ²ₜ
    /// </summary>
    public class GarchModel
    {
        #region 数学常量定义
        
        /// <summary>常数项初值</summary>
        public const double GARCH_OMEGA_DEFAULT = 0.000001;
        
        /// <summary>ARCH项系数初值</summary>
        public const double GARCH_ALPHA_DEFAULT = 0.1;
        
        /// <summary>GARCH项系数初值</summary>
        public const double GARCH_BETA_DEFAULT = 0.85;
        
        /// <summary>波动率阈值系数</summary>
        public const double VOLATILITY_THRESHOLD_FACTOR = 0.8;
        
        /// <summary>最大似然估计迭代次数</summary>
        public const int MAX_LIKELIHOOD_ITERATIONS = 1000;
        
        /// <summary>收敛判断阈值</summary>
        public const double CONVERGENCE_TOLERANCE = 1e-8;
        
        /// <summary>最小样本要求（约1年交易日）</summary>
        public const int MIN_SAMPLE_SIZE = 252;

        /// <summary>最小数据点数量</summary>
        public const int MIN_DATA_POINTS = 20;

        /// <summary>预测时间范围（天）</summary>
        public const int PREDICTION_HORIZON = 5;
        
        /// <summary>计算精度常量</summary>
        public const double CALCULATION_PRECISION = 1e-12;
        
        #endregion

        #region 模型参数和状态
        
        /// <summary>常数项参数 ω</summary>
        public double Omega { get; private set; }
        
        /// <summary>ARCH项系数 α</summary>
        public double Alpha { get; private set; }
        
        /// <summary>GARCH项系数 β</summary>
        public double Beta { get; private set; }
        
        /// <summary>均值参数 μ</summary>
        public double Mu { get; private set; }
        
        /// <summary>长期方差</summary>
        public double LongTermVariance => Omega / (1.0 - Alpha - Beta);
        
        /// <summary>是否已完成参数估计</summary>
        public bool IsEstimated { get; private set; }
        
        /// <summary>模型对数似然值</summary>
        public double LogLikelihood { get; private set; }
        
        /// <summary>收敛迭代次数</summary>
        public int ConvergenceIterations { get; private set; }
        
        #endregion

        #region 构造函数
        
        /// <summary>
        /// 初始化GARCH模型，使用默认参数值
        /// </summary>
        public GarchModel()
        {
            Omega = GARCH_OMEGA_DEFAULT;
            Alpha = GARCH_ALPHA_DEFAULT;
            Beta = GARCH_BETA_DEFAULT;
            Mu = 0.0;
            IsEstimated = false;
            LogLikelihood = double.NaN;
            ConvergenceIterations = 0;
        }
        
        /// <summary>
        /// 使用指定初始参数初始化GARCH模型
        /// </summary>
        /// <param name="omega">常数项初值</param>
        /// <param name="alpha">ARCH项系数初值</param>
        /// <param name="beta">GARCH项系数初值</param>
        /// <param name="mu">均值参数初值</param>
        public GarchModel(double omega, double alpha, double beta, double mu = 0.0)
        {
            if (!ValidateParameters(omega, alpha, beta))
            {
                throw new ArgumentException("GARCH参数不满足约束条件：ω > 0, α ≥ 0, β ≥ 0, α + β < 1");
            }
            
            Omega = omega;
            Alpha = alpha;
            Beta = beta;
            Mu = mu;
            IsEstimated = false;
            LogLikelihood = double.NaN;
            ConvergenceIterations = 0;
        }
        
        #endregion

        #region 参数估计主要方法
        
        /// <summary>
        /// 使用最大似然估计方法估计GARCH模型参数
        /// LL(θ) = -1/2 × Σ[log(2π) + log(σ²ₜ) + εₜ²/σ²ₜ]
        /// </summary>
        /// <param name="returns">收益率时间序列</param>
        /// <returns>参数估计是否成功</returns>
        public bool EstimateParameters(double[] returns)
        {
            if (!ValidateInputData(returns))
            {
                return false;
            }

            try
            {
                var parameters = MaximumLikelihoodEstimation(returns);
                
                if (parameters != null && ValidateParameters(parameters[0], parameters[1], parameters[2]))
                {
                    Omega = parameters[0];
                    Alpha = parameters[1];
                    Beta = parameters[2];
                    Mu = parameters[3];
                    
                    LogLikelihood = CalculateLogLikelihood(returns, parameters);
                    IsEstimated = true;
                    
                    return true;
                }
            }
            catch (Exception ex)
            {
                // 记录异常信息但不抛出，返回估计失败
                Console.WriteLine($"GARCH参数估计异常: {ex.Message}");
            }
            
            return false;
        }
        
        /// <summary>
        /// 最大似然估计核心算法实现
        /// 使用BFGS优化算法进行参数估计
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <returns>估计得到的参数向量 [ω, α, β, μ]</returns>
        private double[] MaximumLikelihoodEstimation(double[] returns)
        {
            // 初始参数设置
            double[] initialParams = GetInitialParameters(returns);
            
            // 参数边界约束
            double[] lowerBounds = { 1e-10, 0.0, 0.0, -1.0 };  // ω > 0, α ≥ 0, β ≥ 0
            double[] upperBounds = { 1.0, 0.999, 0.999, 1.0 }; // 合理上界
            
            // 创建BFGS优化器
            var optimizer = new BfgsOptimizer();
            
            // 定义目标函数（负对数似然，因为我们要最小化）
            Func<double[], double> objectiveFunc = (params_) =>
            {
                // 检查参数约束
                if (!ValidateParameters(params_[0], params_[1], params_[2]))
                {
                    return 1e10; // 返回大的惩罚值
                }
                
                // 返回负对数似然（最小化负值等同于最大化原值）
                return -CalculateLogLikelihood(returns, params_);
            };
            
            // 执行优化
            var result = optimizer.Minimize(
                objectiveFunc,
                initialParams,
                lowerBounds,
                upperBounds,
                MAX_LIKELIHOOD_ITERATIONS,
                CONVERGENCE_TOLERANCE
            );
            
            ConvergenceIterations = result.Iterations;
            
            // 验证结果并返回
            if (result.Converged && ValidateParameters(result.Parameters[0], result.Parameters[1], result.Parameters[2]))
            {
                return result.Parameters;
            }
            
            // 如果BFGS失败，尝试备用方法
            return FallbackEstimation(returns, initialParams);
        }
        
        /// <summary>
        /// 获取初始参数估计
        /// 基于矩估计法提供合理的初始值
        /// </summary>
        private double[] GetInitialParameters(double[] returns)
        {
            // 计算样本统计量
            double mean = returns.Average();
            double variance = CalculateSampleVariance(returns);
            
            // 使用矩估计法设置初始参数
            double omega = variance * 0.05;  // 长期方差的一部分
            double alpha = 0.05;              // 经验值
            double beta = 0.90;               // 高持续性是常见的
            
            // 确保满足稳定性条件
            if (alpha + beta >= 1.0)
            {
                double scale = 0.95 / (alpha + beta);
                alpha *= scale;
                beta *= scale;
            }
            
            return new double[] { omega, alpha, beta, mean };
        }
        
        /// <summary>
        /// 备用估计方法（当BFGS失败时）
        /// 使用网格搜索结合局部优化
        /// </summary>
        private double[] FallbackEstimation(double[] returns, double[] initialParams)
        {
            double[] bestParams = (double[])initialParams.Clone();
            double bestLL = double.NegativeInfinity;
            
            // 网格搜索
            double[] alphaGrid = { 0.01, 0.05, 0.1, 0.15 };
            double[] betaGrid = { 0.7, 0.8, 0.85, 0.9 };
            
            foreach (double alpha in alphaGrid)
            {
                foreach (double beta in betaGrid)
                {
                    if (alpha + beta >= 0.99) continue;
                    
                    double omega = CalculateSampleVariance(returns) * (1 - alpha - beta);
                    double mu = returns.Average();
                    double[] testParams = { omega, alpha, beta, mu };
                    
                    double ll = CalculateLogLikelihood(returns, testParams);
                    if (ll > bestLL)
                    {
                        bestLL = ll;
                        Array.Copy(testParams, bestParams, testParams.Length);
                    }
                }
            }
            
            // 局部优化（简单梯度下降）
            for (int iter = 0; iter < 100; iter++)
            {
                if (!UpdateParameters(bestParams, returns, iter))
                    break;
            }
            
            return ValidateParameters(bestParams[0], bestParams[1], bestParams[2]) 
                ? bestParams : null;
        }
        
        /// <summary>
        /// 参数更新策略（简化版梯度下降）
        /// </summary>
        /// <param name="parameters">当前参数</param>
        /// <param name="returns">收益率序列</param>
        /// <param name="iteration">迭代次数</param>
        /// <returns>是否继续迭代</returns>
        private bool UpdateParameters(double[] parameters, double[] returns, int iteration)
        {
            double learningRate = 0.001 / (1.0 + iteration * 0.01);
            double[] gradient = ComputeGradient(returns, parameters);
            
            bool hasChanged = false;
            for (int i = 0; i < parameters.Length; i++)
            {
                double oldValue = parameters[i];
                parameters[i] += learningRate * gradient[i];
                
                // 应用参数约束
                ApplyParameterConstraints(parameters);
                
                if (Math.Abs(parameters[i] - oldValue) > CONVERGENCE_TOLERANCE)
                {
                    hasChanged = true;
                }
            }
            
            return hasChanged;
        }
        
        /// <summary>
        /// 计算对数似然函数的数值梯度
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <param name="parameters">当前参数</param>
        /// <returns>梯度向量</returns>
        private double[] ComputeGradient(double[] returns, double[] parameters)
        {
            double[] gradient = new double[4];
            double epsilon = 1e-6;
            
            for (int i = 0; i < parameters.Length; i++)
            {
                double[] paramsPlus = (double[])parameters.Clone();
                double[] paramsMinus = (double[])parameters.Clone();
                
                paramsPlus[i] += epsilon;
                paramsMinus[i] -= epsilon;
                
                double llPlus = CalculateLogLikelihood(returns, paramsPlus);
                double llMinus = CalculateLogLikelihood(returns, paramsMinus);
                
                gradient[i] = (llPlus - llMinus) / (2.0 * epsilon);
            }
            
            return gradient;
        }
        
        #endregion

        #region 波动率计算和预测
        
        /// <summary>
        /// 计算对数似然函数值
        /// LL(θ) = -1/2 × Σ[log(2π) + log(σ²ₜ) + εₜ²/σ²ₜ]
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <param name="parameters">参数向量 [ω, α, β, μ]</param>
        /// <returns>对数似然值</returns>
        public double CalculateLogLikelihood(double[] returns, double[] parameters)
        {
            double omega = parameters[0];
            double alpha = parameters[1];
            double beta = parameters[2];
            double mu = parameters[3];
            
            double[] variances = CalculateConditionalVariances(returns, omega, alpha, beta, mu);
            double sumLogLikelihood = 0.0;
            
            for (int t = 1; t < returns.Length; t++)
            {
                double epsilon = returns[t] - mu;
                double variance = Math.Max(variances[t], CALCULATION_PRECISION);
                
                // 累加每个观测的对数似然贡献
                sumLogLikelihood += Math.Log(2.0 * Math.PI) + Math.Log(variance) + 
                                   (epsilon * epsilon) / variance;
            }
            
            // 最终应用 -1/2 系数
            return -0.5 * sumLogLikelihood;
        }
        
        /// <summary>
        /// 计算条件方差序列
        /// σ²ₜ = ω + α × εₜ₋₁² + β × σ²ₜ₋₁
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <param name="omega">常数项</param>
        /// <param name="alpha">ARCH系数</param>
        /// <param name="beta">GARCH系数</param>
        /// <param name="mu">均值</param>
        /// <returns>条件方差序列</returns>
        public double[] CalculateConditionalVariances(double[] returns, double omega, 
            double alpha, double beta, double mu)
        {
            int n = returns.Length;
            double[] variances = new double[n];
            
            // 初始方差设为样本方差
            variances[0] = CalculateSampleVariance(returns);
            
            for (int t = 1; t < n; t++)
            {
                double epsilon_prev = returns[t-1] - mu;
                variances[t] = omega + alpha * epsilon_prev * epsilon_prev + 
                              beta * variances[t-1];
                
                // 确保方差为正值
                variances[t] = Math.Max(variances[t], CALCULATION_PRECISION);
            }
            
            return variances;
        }
        
        /// <summary>
        /// 波动率多步预测
        /// σ²ₜ₊ₕ = σ²_LR + (α+β)ʰ × (σ²ₜ - σ²_LR)
        /// </summary>
        /// <param name="currentVariance">当前方差</param>
        /// <param name="horizons">预测步数数组</param>
        /// <returns>预测方差序列</returns>
        public double[] PredictVolatility(double currentVariance, int[] horizons)
        {
            if (!IsEstimated)
            {
                throw new InvalidOperationException("模型参数尚未估计，无法进行预测");
            }
            
            double longTermVar = LongTermVariance;
            double persistence = Alpha + Beta;
            double[] predictions = new double[horizons.Length];
            
            for (int i = 0; i < horizons.Length; i++)
            {
                int h = horizons[i];
                predictions[i] = longTermVar + Math.Pow(persistence, h) * 
                                (currentVariance - longTermVar);
            }
            
            return predictions;
        }
        
        /// <summary>
        /// 单步波动率预测
        /// </summary>
        /// <param name="currentVariance">当前方差</param>
        /// <returns>下一期预测方差</returns>
        public double PredictNextVolatility(double currentVariance)
        {
            if (!IsEstimated)
            {
                throw new InvalidOperationException("模型参数尚未估计，无法进行预测");
            }

            return Omega + Alpha * currentVariance + Beta * currentVariance;
        }

        /// <summary>
        /// 获取当前波动率
        /// 基于历史数据计算当前条件波动率
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <returns>当前波动率（标准差形式）</returns>
        public double GetCurrentVolatility(double[] returns)
        {
            if (returns == null || returns.Length < MIN_DATA_POINTS)
            {
                throw new ArgumentException($"计算当前波动率需要至少{MIN_DATA_POINTS}个数据点");
            }

            // 如果模型已估计，使用条件方差
            if (IsEstimated)
            {
                var conditionalVariances = CalculateConditionalVariances(
                    returns, Omega, Alpha, Beta, Mu);
                double currentVariance = conditionalVariances[conditionalVariances.Length - 1];
                return Math.Sqrt(currentVariance);
            }
            else
            {
                // 否则使用历史波动率
                double mean = returns.Average();
                double sumSquares = returns.Sum(r => Math.Pow(r - mean, 2));
                double variance = sumSquares / (returns.Length - 1);
                return Math.Sqrt(variance);
            }
        }

        #endregion

        #region 波动率聚类检测
        
        /// <summary>
        /// 检测波动率聚类
        /// 使用K-means聚类和统计分析识别波动率聚类模式
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <param name="windowSize">滚动窗口大小（默认20）</param>
        /// <returns>波动率聚类分析结果</returns>
        public VolatilityCluster DetectVolatilityClustering(double[] returns, int windowSize = 20)
        {
            if (!IsEstimated)
            {
                throw new InvalidOperationException("模型参数尚未估计，无法检测波动率聚类");
            }
            
            // 计算滚动窗口波动率
            double[] rollingVolatility = CalculateRollingVolatility(returns, windowSize);
            
            // 使用K-means聚类（低、中、高三个聚类）
            var kmeans = new KMeansClustering(3);
            int[] clusterLabels = kmeans.Cluster(rollingVolatility);
            double[] centroids = kmeans.GetCentroids();
            
            // 排序中心点以确定低、中、高分类
            var sortedCentroids = centroids.OrderBy(c => c).ToArray();
            int lowLabel = Array.IndexOf(centroids, sortedCentroids[0]);
            int mediumLabel = Array.IndexOf(centroids, sortedCentroids[1]);
            int highLabel = Array.IndexOf(centroids, sortedCentroids[2]);
            
            // 创建聚类结果
            var cluster = new VolatilityCluster();
            
            // 识别各波动率时期
            IdentifyVolatilityPeriods(rollingVolatility, clusterLabels, 
                lowLabel, mediumLabel, highLabel, cluster);
            
            // 识别状态转换点
            IdentifyRegimeChangePoints(clusterLabels, cluster);
            
            // 计算聚类强度
            cluster.ClusteringStrength = CalculateClusteringStrength(
                rollingVolatility, clusterLabels, centroids);
            
            // 计算波动率持续性（自相关）
            cluster.VolatilityPersistence = CalculateVolatilityPersistence(rollingVolatility);
            
            // 计算平均聚类持续时间
            cluster.AverageClusterDuration = CalculateAverageClusterDuration(cluster);
            
            // 确定当前状态
            if (clusterLabels.Length > 0)
            {
                int lastLabel = clusterLabels[clusterLabels.Length - 1];
                cluster.CurrentState = DetermineVolatilityState(lastLabel, lowLabel, mediumLabel, highLabel);
            }
            
            return cluster;
        }
        
        /// <summary>
        /// 计算滚动窗口波动率
        /// </summary>
        private double[] CalculateRollingVolatility(double[] returns, int windowSize)
        {
            int n = returns.Length - windowSize + 1;
            if (n <= 0)
            {
                return new double[] { CalculateSampleVariance(returns) };
            }
            
            double[] volatility = new double[n];
            
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                double sumSq = 0;
                
                for (int j = 0; j < windowSize; j++)
                {
                    double value = returns[i + j];
                    sum += value;
                    sumSq += value * value;
                }
                
                double mean = sum / windowSize;
                double variance = (sumSq / windowSize) - (mean * mean);
                volatility[i] = Math.Sqrt(Math.Max(variance, 0));
            }
            
            return volatility;
        }
        
        /// <summary>
        /// 识别波动率时期
        /// </summary>
        private void IdentifyVolatilityPeriods(
            double[] volatility,
            int[] labels,
            int lowLabel,
            int mediumLabel,
            int highLabel,
            VolatilityCluster cluster)
        {
            int currentLabel = -1;
            int startIndex = 0;
            
            for (int i = 0; i <= labels.Length; i++)
            {
                bool isEnd = (i == labels.Length);
                int label = isEnd ? -1 : labels[i];
                
                if (label != currentLabel)
                {
                    if (currentLabel >= 0 && i > startIndex)
                    {
                        // 创建波动率时期
                        var period = CreateVolatilityPeriod(
                            volatility, startIndex, i - 1);
                        
                        // 添加到相应列表
                        if (currentLabel == lowLabel)
                            cluster.LowVolatilityPeriods.Add(period);
                        else if (currentLabel == mediumLabel)
                            cluster.MediumVolatilityPeriods.Add(period);
                        else if (currentLabel == highLabel)
                            cluster.HighVolatilityPeriods.Add(period);
                    }
                    
                    currentLabel = label;
                    startIndex = i;
                }
            }
        }
        
        /// <summary>
        /// 创建波动率时期对象
        /// </summary>
        private VolatilityPeriod CreateVolatilityPeriod(
            double[] volatility,
            int startIndex,
            int endIndex)
        {
            var period = new VolatilityPeriod
            {
                StartIndex = startIndex,
                EndIndex = endIndex
            };
            
            double sum = 0;
            double min = double.MaxValue;
            double max = double.MinValue;
            int count = 0;
            
            for (int i = startIndex; i <= endIndex && i < volatility.Length; i++)
            {
                double v = volatility[i];
                sum += v;
                min = Math.Min(min, v);
                max = Math.Max(max, v);
                count++;
            }
            
            period.AverageVolatility = count > 0 ? sum / count : 0;
            period.MinVolatility = min;
            period.MaxVolatility = max;
            
            // 计算标准差
            double sumSq = 0;
            for (int i = startIndex; i <= endIndex && i < volatility.Length; i++)
            {
                double diff = volatility[i] - period.AverageVolatility;
                sumSq += diff * diff;
            }
            period.VolatilityStdDev = count > 1 ? 
                Math.Sqrt(sumSq / (count - 1)) : 0;
            
            return period;
        }
        
        /// <summary>
        /// 识别状态转换点
        /// </summary>
        private void IdentifyRegimeChangePoints(int[] labels, VolatilityCluster cluster)
        {
            for (int i = 1; i < labels.Length; i++)
            {
                if (labels[i] != labels[i - 1])
                {
                    cluster.RegimeChangePoints.Add(i);
                }
            }
        }
        
        /// <summary>
        /// 计算聚类强度
        /// 使用轮廓系数（Silhouette Score）衡量聚类质量
        /// </summary>
        private double CalculateClusteringStrength(
            double[] volatility,
            int[] labels,
            double[] centroids)
        {
            if (volatility.Length <= 1 || centroids.Length <= 1)
                return 0;
            
            double totalScore = 0;
            int count = 0;
            
            for (int i = 0; i < volatility.Length; i++)
            {
                // 计算到自己聚类中心的平均距离
                double a = Math.Abs(volatility[i] - centroids[labels[i]]);
                
                // 计算到最近的其他聚类中心的平均距离
                double b = double.MaxValue;
                for (int k = 0; k < centroids.Length; k++)
                {
                    if (k != labels[i])
                    {
                        double dist = Math.Abs(volatility[i] - centroids[k]);
                        b = Math.Min(b, dist);
                    }
                }
                
                // 计算轮廓系数
                double s = (b - a) / Math.Max(a, b);
                totalScore += s;
                count++;
            }
            
            // 归一化到[0,1]
            double avgScore = count > 0 ? totalScore / count : 0;
            return Math.Max(0, Math.Min(1, (avgScore + 1) / 2));
        }
        
        /// <summary>
        /// 计算波动率持续性（一阶自相关）
        /// </summary>
        private double CalculateVolatilityPersistence(double[] volatility)
        {
            if (volatility.Length < 2)
                return 0;
            
            double mean = volatility.Average();
            double cov = 0;
            double var = 0;
            
            for (int i = 1; i < volatility.Length; i++)
            {
                cov += (volatility[i] - mean) * (volatility[i - 1] - mean);
            }
            
            for (int i = 0; i < volatility.Length; i++)
            {
                var += (volatility[i] - mean) * (volatility[i] - mean);
            }
            
            return var > 0 ? cov / var : 0;
        }
        
        /// <summary>
        /// 计算平均聚类持续时间
        /// </summary>
        private double CalculateAverageClusterDuration(VolatilityCluster cluster)
        {
            var allPeriods = new List<VolatilityPeriod>();
            allPeriods.AddRange(cluster.LowVolatilityPeriods);
            allPeriods.AddRange(cluster.MediumVolatilityPeriods);
            allPeriods.AddRange(cluster.HighVolatilityPeriods);
            
            if (allPeriods.Count == 0)
                return 0;
            
            return allPeriods.Average(p => p.Duration);
        }
        
        /// <summary>
        /// 确定波动率状态
        /// </summary>
        private VolatilityState DetermineVolatilityState(
            int label, int lowLabel, int mediumLabel, int highLabel)
        {
            if (label == lowLabel)
                return VolatilityState.Low;
            else if (label == mediumLabel)
                return VolatilityState.Medium;
            else if (label == highLabel)
                return VolatilityState.High;
            else
                return VolatilityState.Transitioning;
        }
        
        #endregion

        #region 模型诊断和验证
        
        /// <summary>
        /// 执行模型诊断检验
        /// </summary>
        /// <param name="returns">原始收益率序列</param>
        /// <returns>诊断结果</returns>
        public GarchDiagnostics PerformDiagnostics(double[] returns)
        {
            if (!IsEstimated)
            {
                throw new InvalidOperationException("模型参数尚未估计，无法执行诊断");
            }
            
            var diagnostics = new GarchDiagnostics();
            
            // 计算标准化残差
            double[] standardizedResiduals = CalculateStandardizedResiduals(returns);
            
            // Ljung-Box检验
            diagnostics.LjungBoxPValue = LjungBoxTest(standardizedResiduals, 10);
            
            // ARCH效应检验
            diagnostics.ArchTestPValue = ArchTest(standardizedResiduals, 5);
            
            // 正态性检验
            diagnostics.NormalityTestPValue = JarqueBeraTest(standardizedResiduals);
            
            // 残差分析
            diagnostics.ResidualMean = standardizedResiduals.Average();
            diagnostics.ResidualStdDev = Math.Sqrt(CalculateSampleVariance(standardizedResiduals));
            
            return diagnostics;
        }
        
        /// <summary>
        /// 计算标准化残差
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <returns>标准化残差序列</returns>
        public double[] CalculateStandardizedResiduals(double[] returns)
        {
            double[] variances = CalculateConditionalVariances(returns, Omega, Alpha, Beta, Mu);
            double[] residuals = new double[returns.Length];
            
            for (int t = 0; t < returns.Length; t++)
            {
                residuals[t] = (returns[t] - Mu) / Math.Sqrt(variances[t]);
            }
            
            return residuals;
        }
        
        /// <summary>
        /// Ljung-Box检验统计量（简化实现）
        /// </summary>
        /// <param name="residuals">残差序列</param>
        /// <param name="lags">滞后阶数</param>
        /// <returns>p值</returns>
        private double LjungBoxTest(double[] residuals, int lags)
        {
            double n = residuals.Length;
            double[] autocorr = CalculateAutocorrelations(residuals, lags);
            double lbStatistic = 0.0;
            
            for (int k = 1; k <= lags; k++)
            {
                lbStatistic += autocorr[k-1] * autocorr[k-1] / (n - k);
            }
            
            lbStatistic *= n * (n + 2);
            
            // 简化的p值计算（实际应使用卡方分布）
            return 1.0 - Math.Min(0.99, lbStatistic / (2.0 * lags));
        }
        
        /// <summary>
        /// ARCH效应检验（简化实现）
        /// </summary>
        /// <param name="residuals">残差序列</param>
        /// <param name="lags">滞后阶数</param>
        /// <returns>p值</returns>
        private double ArchTest(double[] residuals, int lags)
        {
            double[] squaredResiduals = residuals.Select(x => x * x).ToArray();
            double[] autocorr = CalculateAutocorrelations(squaredResiduals, lags);
            
            double archStatistic = 0.0;
            for (int k = 0; k < lags; k++)
            {
                archStatistic += autocorr[k] * autocorr[k];
            }
            
            archStatistic *= residuals.Length;
            
            // 简化的p值计算
            return 1.0 - Math.Min(0.99, archStatistic / (2.0 * lags));
        }
        
        /// <summary>
        /// Jarque-Bera正态性检验（简化实现）
        /// </summary>
        /// <param name="residuals">残差序列</param>
        /// <returns>p值</returns>
        private double JarqueBeraTest(double[] residuals)
        {
            double n = residuals.Length;
            double mean = residuals.Average();
            double variance = CalculateSampleVariance(residuals);
            double stdDev = Math.Sqrt(variance);
            
            // 计算偏度和峰度
            double skewness = 0.0, kurtosis = 0.0;
            foreach (double x in residuals)
            {
                double standardized = (x - mean) / stdDev;
                double pow2 = standardized * standardized;
                double pow3 = pow2 * standardized;
                double pow4 = pow2 * pow2;
                
                skewness += pow3;
                kurtosis += pow4;
            }
            
            skewness /= n;
            kurtosis = kurtosis / n - 3.0; // 超额峰度
            
            double jbStatistic = (n / 6.0) * (skewness * skewness + 0.25 * kurtosis * kurtosis);
            
            // 简化的p值计算
            return 1.0 - Math.Min(0.99, jbStatistic / 10.0);
        }
        
        #endregion

        #region 辅助方法
        
        /// <summary>
        /// 验证模型参数约束条件
        /// ω > 0, α ≥ 0, β ≥ 0, α + β < 1
        /// </summary>
        /// <param name="omega">常数项</param>
        /// <param name="alpha">ARCH系数</param>
        /// <param name="beta">GARCH系数</param>
        /// <returns>参数是否满足约束</returns>
        public static bool ValidateParameters(double omega, double alpha, double beta)
        {
            return omega > 0 && alpha >= 0 && beta >= 0 && (alpha + beta) < 1.0;
        }
        
        /// <summary>
        /// 应用参数约束
        /// </summary>
        /// <param name="parameters">参数数组</param>
        private void ApplyParameterConstraints(double[] parameters)
        {
            // 确保 ω > 0
            parameters[0] = Math.Max(parameters[0], CALCULATION_PRECISION);
            
            // 确保 α ≥ 0
            parameters[1] = Math.Max(parameters[1], 0.0);
            
            // 确保 β ≥ 0
            parameters[2] = Math.Max(parameters[2], 0.0);
            
            // 确保 α + β < 1
            if (parameters[1] + parameters[2] >= 1.0)
            {
                double sum = parameters[1] + parameters[2];
                parameters[1] = parameters[1] / sum * 0.99;
                parameters[2] = parameters[2] / sum * 0.99;
            }
        }
        
        /// <summary>
        /// 验证输入数据有效性
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <returns>数据是否有效</returns>
        private bool ValidateInputData(double[] returns)
        {
            return returns != null && returns.Length >= MIN_SAMPLE_SIZE && 
                   returns.All(x => !double.IsNaN(x) && !double.IsInfinity(x));
        }
        
        /// <summary>
        /// 计算样本方差
        /// </summary>
        /// <param name="data">数据序列</param>
        /// <returns>样本方差</returns>
        private double CalculateSampleVariance(double[] data)
        {
            double mean = data.Average();
            double variance = data.Select(x => (x - mean) * (x - mean)).Average();
            return variance;
        }
        
        /// <summary>
        /// 计算自相关函数
        /// </summary>
        /// <param name="data">数据序列</param>
        /// <param name="maxLag">最大滞后阶数</param>
        /// <returns>自相关系数数组</returns>
        private double[] CalculateAutocorrelations(double[] data, int maxLag)
        {
            int n = data.Length;
            double mean = data.Average();
            double[] autocorr = new double[maxLag];
            double c0 = 0.0;
            
            // 计算0阶自协方差
            for (int t = 0; t < n; t++)
            {
                c0 += (data[t] - mean) * (data[t] - mean);
            }
            c0 /= n;
            
            // 计算各阶自相关系数
            for (int k = 1; k <= maxLag; k++)
            {
                double ck = 0.0;
                for (int t = k; t < n; t++)
                {
                    ck += (data[t] - mean) * (data[t-k] - mean);
                }
                ck /= n;
                autocorr[k-1] = ck / c0;
            }
            
            return autocorr;
        }
        
        /// <summary>
        /// 获取模型摘要信息
        /// </summary>
        /// <returns>模型摘要字符串</returns>
        public string GetModelSummary()
        {
            if (!IsEstimated)
            {
                return "GARCH(1,1)模型 - 参数尚未估计";
            }
            
            return $"GARCH(1,1)模型摘要:\n" +
                   $"ω (常数项) = {Omega:F8}\n" +
                   $"α (ARCH系数) = {Alpha:F6}\n" +
                   $"β (GARCH系数) = {Beta:F6}\n" +
                   $"μ (均值) = {Mu:F6}\n" +
                   $"长期方差 = {LongTermVariance:F8}\n" +
                   $"持续性 (α+β) = {Alpha + Beta:F6}\n" +
                   $"对数似然值 = {LogLikelihood:F4}\n" +
                   $"收敛迭代次数 = {ConvergenceIterations}";
        }
        
        #endregion

        #region 辅助方法和兼容性接口

        /// <summary>
        /// 获取估计的参数数组
        /// </summary>
        /// <returns>参数数组 [omega, alpha, beta]</returns>
        public double[] GetEstimatedParameters()
        {
            return new double[] { Omega, Alpha, Beta };
        }

        /// <summary>
        /// 拟合模型（EstimateParameters的别名）
        /// </summary>
        /// <param name="returns">收益率序列</param>
        /// <returns>是否拟合成功</returns>
        public bool FitModel(double[] returns)
        {
            return EstimateParameters(returns);
        }

        /// <summary>
        /// Fit方法（兼容性）
        /// </summary>
        public bool Fit(double[] returns)
        {
            return EstimateParameters(returns);
        }

        /// <summary>
        /// 获取无条件方差
        /// </summary>
        /// <returns>无条件方差值</returns>
        public double GetUnconditionalVariance()
        {
            if (!IsEstimated)
            {
                return double.NaN;
            }

            // 无条件方差 = ω / (1 - α - β)
            double denominator = 1.0 - Alpha - Beta;
            if (denominator <= 0)
            {
                return double.NaN;
            }

            return Omega / denominator;
        }

        /// <summary>
        /// UnconditionalVariance属性（兼容性）
        /// </summary>
        public double UnconditionalVariance => GetUnconditionalVariance();

        /// <summary>
        /// 重载的PredictVolatility方法（单期预测）
        /// </summary>
        /// <param name="horizons">预测期数</param>
        /// <returns>预测的波动率数组</returns>
        public double[] PredictVolatility(int horizons)
        {
            if (!IsEstimated)
            {
                throw new InvalidOperationException("模型尚未估计，无法预测");
            }

            // 使用最后的条件方差作为起始点
            double lastVariance = GetUnconditionalVariance();
            return PredictVolatility(Math.Sqrt(lastVariance), new int[] { horizons });
        }

        #endregion
    }

    #region 辅助数据结构
    
    /// <summary>
    /// GARCH模型诊断结果
    /// </summary>
    public class GarchDiagnostics
    {
        /// <summary>Ljung-Box检验p值</summary>
        public double LjungBoxPValue { get; set; }
        
        /// <summary>ARCH效应检验p值</summary>
        public double ArchTestPValue { get; set; }
        
        /// <summary>正态性检验p值</summary>
        public double NormalityTestPValue { get; set; }
        
        /// <summary>残差均值</summary>
        public double ResidualMean { get; set; }
        
        /// <summary>残差标准差</summary>
        public double ResidualStdDev { get; set; }
        
        /// <summary>
        /// 检查模型是否通过诊断检验
        /// </summary>
        /// <param name="significanceLevel">显著性水平</param>
        /// <returns>是否通过检验</returns>
        public bool PassesDiagnostics(double significanceLevel = 0.05)
        {
            return LjungBoxPValue > significanceLevel && 
                   ArchTestPValue > significanceLevel &&
                   Math.Abs(ResidualMean) < 0.1 &&
                   Math.Abs(ResidualStdDev - 1.0) < 0.2;
        }
        
        /// <summary>
        /// 获取诊断摘要
        /// </summary>
        /// <returns>诊断摘要字符串</returns>
        public string GetDiagnosticSummary()
        {
            return $"GARCH模型诊断结果:\n" +
                   $"Ljung-Box检验 p值 = {LjungBoxPValue:F4} {(LjungBoxPValue > 0.05 ? "✓" : "✗")}\n" +
                   $"ARCH效应检验 p值 = {ArchTestPValue:F4} {(ArchTestPValue > 0.05 ? "✓" : "✗")}\n" +
                   $"正态性检验 p值 = {NormalityTestPValue:F4} {(NormalityTestPValue > 0.05 ? "✓" : "✗")}\n" +
                   $"残差均值 = {ResidualMean:F6}\n" +
                   $"残差标准差 = {ResidualStdDev:F6}\n" +
                   $"整体诊断: {(PassesDiagnostics() ? "通过 ✓" : "未通过 ✗")}";
        }

        #endregion
    }
}
